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This  is  the  first  report  of  the  project  to  apply  geostatistics  and  wavelet  analysis  toj 
Dart  of  the  area  chosen  at  Fort  A.P.  Hill  in  northeastern  Virginia.  This  reportj 
smbraces  a  summary  of  time  spent  at  TEC,  cokriging  of  Korean  temperature! 
data, -and  the  comparison  between  geostatistics  and  wavelet  analysis  in  these 
sections.  Although  this  report  is  the  first,  it  covers  almost  half  of  the  work  to  be 
done.  The  cokriging  of  temperature  in  Korea  is  an  exercise  to  determine  whether 
sstimates  can  be  improved  by  using  more  information  on  altitude  to  estimate 
:emperature  with  smaller  error  than  by  ordinary  kriging.  Wavelets  and  geostatistics 
:an  both  be  used  for  filtering  data  for  data  reconstruction.  A  comparison  between 
:he  two  approaches  is  described.  ■  It  seems  that  wavelets  provide  a  more  accurate 
nethod  for  data  reconstruction,  but  geostatistics  is  more  appropriate  for  exploring 
different  resolutions  of  spatial  variation  that  have  been  identified  by  the  variogram. 
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Evaluation  of  geostatistics  and  wavelets  for  identifying  relations  between  imagery 
and  different  spatial  resolutions  and  for  data  compression 


Introduction 

This  report  embraces  several  aspects  of  the  work  to  date.  It  begins  with  a  brief 
summary  of  the  work  that  was  done  by  Dr  Oliver  at  TEC  in  February  1997  which  was 
part  of  this  contract  (albeit  slightly  premature),  and  the  work  done  while  E.  Bosch  was 
working  with  Dr  Oliver  at  Reading  in  September.  The  majority  of  the  report  will  focus 
on  a  comparison  between  wavelet  and  geostatistics  as  far  as  we  have  gone. 

The  work  for  most  of  this  project  will  be  based  at  Fort  A.  P.  Hill  in  northeastern 
Virginia,  about  75  miles  from  Washington,  DC.  The  area  is  intensely  dissected  by 
many  small  waterways,  and  this  appears  to  have  contributed  to  the  pattern  of  variation 
observed  in  the  image.  The  ground  cover  survey  data  and  digital  elevation  model  will 
be  used  in  a  second  stage  of  analysis. 


Report  of  visit  to  TEC  in  January  1998 

Much  of  the  first  day  at  TEC  was  spent  discussing  the  results  of  the  first  analyses  from 
Fort  A.  P.  Hill,  and  what  other  work  should  be  done.  In  addition  the  paper  that  has 
now  been  accepted  by  the  International  Journal  of  Remote  Sensing  was  also  discussed 
and  suggestions  for  improvement  considered  and  incorporated.  Since  Dr  Oliver  was  to 
brief  the  senior  management  team  at  TEC  including  Dr  Roper  the  contents  of  the 
briefing  were  also  ratified  at  the  outset. 

The  main  aim  of  this  visit  was  to  work  with  Jim  Shine  to  enable  him  to  make  full  use  of 
Genstat.  A  set  of  programs  was  prepared  to  cover  exploratory  data  analysis 
(histograms,  box  plots,  summary  statistics,  trend  detection  and  so  on),  variogram 
analysis  and  modelling  and  kriging.  All  of  the  programs  were  examined.  They  were 
eventually  compiled  as  part  of  the  aide  memoire  that  formed  an  Aapendix  of  the  final 
report  for  the  previous  contract  (Contract  N68171-97-C-9029). 

Jim  Shine  and  Dr  Oliver  worked  though  all  of  the  programs.  A  problem  was  identified 
with  the  kriging  algorithm  in  Genstat  which  was  eventually  reported  to  the  NAG 
library  and  corrected.  TEC  then  received  a  new  implementation  of  the  package.  At 
least  half  of  the  time  at  TEC  was  spent  instructing  Jim  in  the  use  of  the  programs  and 
interpreting  the  results.  In  addition  we  had  several  discussions  on  geostatistics. 

Half  a  day  was  spent  on  the  briefing  to  Dr  Roper  and  senior  staff  at  TEC,  and  in 
answering  questions  arising  from  this.  During  the  course  of  our  collaboration  we  have 
covered  a  substantial  amount  of  work  and  much  of  it  was  described  briefly  at  this 
meeting.  Dr  Roper  showed  considerable  interest  in  what  has  been  done  and  when  he 
visited  the  University  of  Reading  in  November  1998  it  was  clear  that  he  had  a  sound 


appreciation  of  the  value  of  geostatistical  analysis.  The  discussion  that  followed  the 
briefing  was  stimulating  and  well  considered. 

Other  discussions  were  held  with  Edward  Bosch  about  comparing  wavelets  and 
geostatistics.  This  culminated  with  an  arrangement  for  him  to  visit  Reading  in 
September  1998. 

Dr  Oliver  and  E.  Bosch  worked  together  for  a  week.  The  time  was  used  for  analyses, 
interpreting  results  and  discussion.  Several  analyses  were  undertaken  -  some  of  which 
feature  in  the  report.  Others  have  been  done  by  both  of  us  subsequently.  The  visit  was 
very  profitable  to  both  of  us.  As  a  result  of  this  investigation  we  have  submitted  an 
abstract  to  the  Geostatistics  Congress  to  be  held  in  April  2000.  This  is  appended  at  the 
end  of  the  report.  The  paper,  if  accepted,  will  acknowledge  the  support  of  US  Army 
and  of  TEC  in  this  work,  and  will  be  authored  jointly  by  M.  A.  Oliver,  E.  Bosch  and  K. 
Slocum. 


Comparing  wavelets  and  kriging  for  exploring  nested  scales  of  variation 

This  work  has  continued  to  use  part  of  the  SPOT  image  around  Anderson  Camp. 
However,  it  is  slightly  different  in  extent  from  that  used  previously  to  accommodate 
the  wavelet  analysis.  This  analysis  works  better  if  the  area  analysed  is  some  factor  of  2; 
the  area  chosen  was  27.  As  a  result  the  variogram  analysis  and  factorial  kriging  had  to 
be  redone  so  that  the  results  relate  to  this  new  region  and  are  comparable  with  all 
subsequent  analyses.  The  theory  of  factorial  kriging  was  given  in  the  final  report  for 
Project  3  (Contract  number  N68171-97-C-9029).  The  theory  of  factorial  kriging  was 
also  given  in  the  paper  presented  at  the  geoENV98  conference  which  has  already  been 
submitted  to  TEC.  I  shall  summarise  the  application  of  wavelets,  but  with  little  theory 
as  this  is  not  my  field  and  Edward  Bosch  is  the  local  expert  at  TEC. 

Summary  of  wavelet  analysis 

The  wavelet  transform  has  some  similarity  with  the  Fourier  and  the  windowed  Fourier 
transforms.  However,  the  wavelet  transform  is  localized  in  terms  of  frequency  and 
scale  (via  dilations)  and  in  time  (via  translations),  whereas  the  Fourier  transform, 
although  localized  in  frequency,  is  not  localized  in  time  or  space.  The  windowed 
Fourier  transform  operates  more  locally.  The  advantage  of  the  wavelet  transform  over 
the  windowed  Fourier  transform  is  that  low  frequency  and  high  frequency  resolutions 
can  be  characterized  simultaneously.  This  means  that  wavelet  analysis  is  suitable  for 
situations  where  there  are  different  levels  of  variation  superimposed  on  each  other 
(Daubechies,  1992).  Wavelets  are  also  good  for  describing  transient  data  whereas  the 
Fourier  transform  is  not.  Wavelet  analysis  is  not  affected  by  local  non-stationarity  and 
this  is  an  advantage  it  has  compared  with  geostatistics,  which  assumes  that  the  data  are 
at  least  quasi-stationary  (i.e.  locally  stationary).  Local  non-stationarity  can  arise  where 
there  are  marked  boundaries  that  result  in  a  marked  change  in  the  local  means  of  the 
variable  of  interest. 


2 


Wavelets  are  oscillatory  components  that  operate  locally.  The  wavelet  analysis  starts 
with  the  choice  of  a  mother  wavelet,  w(t),  which  is  fixed.  The  mother  wavelet  can  be 
dilated  or  shrunk  to  examine  components  in  the  variation  that  occur  at  different  spatial 
or  temporal  scales.  This  enables  multi-resolution  analysis  where  different  levels  of 
variation  are  superimposed  on  one  another  (Mallat,  1998).  This  is  our  first  aim  in  this 
investigation.  Redundancy  is  a  major  problem  with  image  data  because  of  the  amount 
of  information  involved.  Wavelets  are  also  of  great  value  for  data  compression  because 
they  are  able  to  remove  and  to  retain  the  important  structure. 


Theory 

Wavelet  analysis  allows  a  signal  (information)  to  be  represented  in  terms  of  a  set  of 
basis  functions,  i.e.  basis  vectors  or  kernels.  The  basis  functions  are  a  set  of  linearly 
independent  functions  that  can  be  used  to  produce  all  admissible  functions  of  /  (t) 
(Strang  and  Nguyen,  1996).  Choosing  the  basis  functions  determines  the  kind  of 
information  that  can  be  extracted. 


fit)  =  Zbsuwsu(t)  (1) 

where  are  the  coefficients  for  scaling  ( s )  and  translation  («).  The  special  feature  of 
the  wavelet  basis  is  that  all  functions  wsu(t)  are  constructed  from  a  single  mother 
wavelet  w(/).  If  the  wavelet  normally  starts  at  t  =  0  and  ends  at  t  —  N,  then  shifted 
wavlets,  wou(t),  start  at  t  =  u  and  end  at  t-  u+N.  The  rescaled  wavelets,  yvscfit),  start  at 
t= 0  and  end  at  t=  N/fi,  and  are  translated  to  the  right  by  s  and  compressed  by  a  factor 
of  '1.  A  typical  wavelet  wsu  is  compressed  j  times  and  shifted  u  times: 

w„(t)  =  w(2‘t-k)  (2) 

Most  wavelets  are  orthogonal  and 

£ Kwsu(t)wsu(t)dt  =  0  (3) 

This  leads  to  a  simple  formula  for  each  coefficient  bSu  in  the  expansion  of  fit).  The 
expansion  in  equation  (1)  is  mutiplied  by  wjk  and  integrated  by: 

£/(0wro  (*)&  =  bsu  £  (wsu  (t)fdt.  (4) 


All  other  terms  disappear  because  of  the  orthogonality. 
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Multiresolution 


The  high  frequency  filter  leads  to  w(/)  and  the  l°w  frequency  one  to  the  scaling 
function  <}>(/).  In  most  wavelet  analyses  the  low  frequency  filter  is  applied  first  and  the 
scaling  function  is  obtained  before  the  wavelet.  The  wavelet  follows  from  <J>(/)  by  one 
application  of  the  high  frequency  filter.  Therefore,  the  scaling  functions  average  or 
smooth  the  data  while  the  detail  in  the  signal  is  retained  in  the  wavelets.  The  scaling 
functions  and  wavelets  comprise  the  components  of  the  multi-resolution  situation 
which  enables  wavelet  analysis  to  examine  features  locally  at  different  scales.  At  a 
given  resolution  the  scaling  functions  4>(2V  -  u)  are  a  basis  for  the  set  of  signals.  The 
level  is  set  by  s,  and  the  steps  at  the  level  are  2's.  The  new  detail  at  the  level  5  is 
represented  by  the  wavelets  w(2*t  -  u). 


Strang  and  Nguyen  (1996)  represent  this  as  follows: 

signal  at  level  s  (local  averages)  N 

+  signal  at  level  s  +1 

details  at  level  s  (local  differences) 

Thus  the  signal,  such  as  the  NIR  information  from  an  image,  can  be  divided  into 
different  scales  of  resolution  by  wavelet  analysis.  Multiresolution  divides  the 
frequencies  into  octave  bands,  from  w  to  2 w,  rather  than  different  frequencies. 

The  construction  of  wavelets  starts  with  vectors.  Every  second  vector  xy  is  a 
combination  of  the  basis  vectors  (1,0)  and  (0,1)  (Strang,  1989).  The  original  data  are 
transformed  by  the  wavelet  coefficients  which  are  both  additive  and  orthogonal.  The 
coefficients  provide  us  with  a  measure  of  the  energy  the  basis  vector  has  at  time  t  or 
scale  s.  The  wavelet  transform  decomposes  the  signal  to  a  basis  (set)  of  elementary 
signals.  It  uses  smaller  windows  at  higher  frequencies  and  larger  ones  at  low 
frequencies.  This  means  that,  in  essence,  the  analysis  is  based  on  a  pair  of  filters  -  one 
of  low  frequency  which  smooths  the  data  and  the  other  is  high  frequency  and  fine 
which  is  related  to  the  wavelet.  The  low  frequency  filter  produces  the  scaling 
coefficients  and  the  high  frequency  one  the  wavelet  coefficients. 


The  wavelet  is  a  function  of  zero  average: 

J"  \|/(r)dt =  0  (^) 

which  is  dilated  or  compressed  with  a  scale  parameter,  s,  and  translated  by  u.  The 
wavelet  coefficient  wflujs)  measures  the  variation  of/ in  a  neighbourhood  u  whose  size 
is  proportional  to  5.  The  method  assumes  a  finite  variance. 

There  are  many  different  kinds  of  wavelets,  the  most  simple  is  the  Haar,  but  probably 
the  most  popular  at  the  moment  are  those  of  Daubechies  (Daubechies,  1988). 
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Analysis  of  the  A.P.  Hill  data 

SPOT  Image 


The  part  of  the  scene  covering  Fort  A.  P.  Hill,  Figure  1,  is  slightly  smaller  than  that 
used  before  (see  report  N68171-97-C-9029),  but  it  covers  the  same  part  of  the  image. 
Analyses  were  carried  out  on  the  complete  data  set  and  on  sub-samples  of  1  pixel  in  2 
for  each  column  and  row  (or  1  pixel  from  a  block  of  4),  1  pixel  in  4  for  each  row  and 
column  (1  pixel  from  a  block  of  16),  and  1  pixel  in  8  for  each  row  and  column  (1  pixel 
from  a  block  of  64  pixels).  The  sub-sets  were  used  to  assess  the  accuracy  of  data 
reconstruction  by  the  two  methods.  Table  1  gives  the  summary  statistics  for  the  full 
data  set. 


Table  1:  Summary  statistics  for  NIR  for  the  128  pixels  by  128  pixels  region  of  Fort  A. 
P.  Hill 


Statistic 

NIR 

Hermite 

polynomials  of  NIR 

Count 

16384.0 

16384.0 

Minimum 

37.00 

-4.496 

Maximum 

183.00 

4.009 

Mean 

117.83 

0.0287 

Variance 

268.99 

0.9995 

Standard  deviation 

16.401 

0.9997 

Skewness 

-0.7408 

-0.0943 

Wavelet  analysis 


The  method  of  wavelet  analysis  that  Edward  Bosch  used  was  that  of  Daubechies 
wavelets  (Daubechies,  1988).  The  size  of  the  image  was  27,  i.e.  128  rows  and  128 
columns  of  pixel  information  for  the  NIR  waveband.  This  region  was  chosen  from  the 
part  of  the  SPOT  image  that  we  analysed  and  described  in  the  previous  report  (1998), 
and  the  size  was  such  to  avoid  any  need  to  pad  the  data  to  create  appropriate 
resolution  levels  for  the  wavelet  analysis. 

The  wavelet  transform  was  done  with  a  pair  of  filters  and  not  the  wavelets  themselves. 
The  signal  is  convolved  by  the  filter  -  the  result  is  that  every  other  point  in  the  output 
of  each  convolution  is  discarded.  The  convolutions  can  be  represented  in  terms  of  two 
matrices  L  and  H.  The  matrix  L  is  made  up  of  shifts  of  the  low  frequency  filter,  and  H 
is  made  up  of  shifts  of  the  high  frequency  filter.  Since  the  wavelets  are  orthogonal: 

H'H  +  L'L  =  I 

where  H'  is  the  transpose  of  H,  L'  is  the  transpose  of  L,  and  I  is  the  identity  matrix. 
The  matrices  H  and  L  correspond  to  the  forward  wavelet,  and  H'  and  1/  to  the  inverse 
wavelet  transform.  If  the  vector  /is  the  frequency  (NIR  in  our  analysis),  then  the  low 
frequency  component  of  the  variation  in  /  is  contained  in  1/  and  the  high  frequency 
component  in  H/. 

L'(  L f)  +  H'(  HO  =  (L'  L V  +(EL'Wjf=  (L'  L)  +  (H'  W  =  1/=/ 

The  vectors  1/  and  H/  contain  half  of  the  number  of  samples  as  in  the  original  set  of 
data  for  /. 

The  filters  were  constructed  so  that  there  were  three  orthogonal  filters  and  three 
vanishing  moments  were  satisfied.  The  more  vanishing  moments  there  are  the  greater 
the  smoothing  (Daubechies,  1992).  With  three  vanishing  moments  the  high  pass  filter 
zeros  out  leaving  behind  little  information.  Most  of  the  ‘energy’  or  information  is  then 
retained  by  the  low  pass  filter.  The  scaling  coefficients  from  the  low  pass  filter  contain 
most  of  the  information  and  it  these  that  were  used  to  reconstruct  the  image. 


If  the  discrete  wavelet  transform  (DWT)  is  applied  once  to  the  data  the  resolution 
becomes  26,  i.e.  resolution  6.  This  results  in  four  quarter  sets  of  data  of  size  64  by  64 
pixels.  The  first  quarter  of  the  data  represent  in  essence  a  sample  of  1  in  2  of  the  rows 
and  columns  of  the  data  matrix.  This  quarter  contains  the  scaling  coefficients  which 
correspond  to  the  low  pass  filter.  The  other  three  quarters  contain  the  wavelet 
coefficients  which  are  high  frequency:  quarter  two  contains  the  vertical  coefficients, 
quarter  3  the  horizontal  ones  and  quarter  4  the  diagonal  coefficients. 
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The  image  information  on  NIR  was  reconstructed  by  an  inverse  wavelet  transform.  For 
this  analysis  the  low  and  high  frequency  components  were  reconstructed  separately. 
For  the  low  frequency  reconstruction: 

JJ=  L'(L/)  +  0'(H/)  =  L'(  L/), 

where  O  is  the  zero  matrix.  Therefore,./?  represents  the  reconstructed  values  of  NIR 
without  the  high  frequency  information.  For  the  high  frequency  reconstruction: 

fh  =  0'(  L/)  +  H'(  H/)  =H'(H/), 


where  fh  is  the  high  frequency  information  only.  Therefore: 


Geostaiistical  analyses 

The  variogram  was  computed  and  modelled  as  usual  using  the  smaller  data  set  with  a 
total  of  16  384  pixels.  Figure  1.  This  variogram  was  then  used  with  the  pixel 
information  to  filter  the  information  by  factorial  kriging  into  the  long-range  and  short- 
range  components.  Ordinary  kriging  was  used  to  estimate  the  values  of  NIR  at 
positions  where  pixels  had  been  removed  from  the  data.  In  other  words  the  estimates 
coincided  with  the  locations  of  the  original  values  so  that  a  direct  comparison  could  be 
made  between  the  estimates  and  these. 


Results 

The  variogram  for  the  new  data  is  still  a  nested  structure,  but  the  correlation  ranges  are 
smaller  than  for  the  larger  part  of  the  scene  that  we  investigated  before.  The  model 
fitted  was  a  nested  spherical  function  with  two  structures.  Since  the  variogram  was 
somewhat  wavy  at  the  longer  lags,  to  improve  the  fit  I  modelled  it  to  a  lag  of  40  only. 
The  short-range  structure  was  6.6  pixels  or  130  m  and  the  long-range  structure  was  21 
pixels  or  420  m.  The  experimental  variogram  (points)  and  the  fitted  model  (line)  are 
given  in  Figure  2  a.  The  parameters  of  the  models  fitted  to  100  and  to  40  lags  are  given 
in  Table  2.  Since  the  data  were  skewed  I  transformed  them  using  Hermite  polynomials 
and  computed  the  variogram  from  the  transformed  data,  Figure  2b  and  Table  2.  There 
was  little  difference  from  the  raw  variogram,  therefore,  I  did  the  analyses  on  the  raw 
data. 
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a) 


Figure  2:  a)  Variogram  of  the  near  infra  read  (NIR)  of  part  of  the  SPOT  image 
(128  by  128  pixels)  for  Fort  A.  P.  Hill,  b)  variogram  of  the  transformed  pixel  data 
using  Hermite  polynomials 
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Multiresolution  analysis  (filtering) 


The  variogram  suggests  that  there  are  two  clear  scales  of  spatial  variation  present:  one 
of  about  120  m  and  the  other  of  about  420.  This  is  also  evident  in  the  pixel  map  of  the 
ordinary  kriged  estimates.  Figure  3.  There  is  local  detailed  variation  superimposed  on  a 
broader  pattern  of  variation.  The  major  large  structures  in  the  variation  that  are  evident 
appear  to  be  related  to  major  relief  forms:  the  drainage  basins  and  the  intervening 
spurs,  and  the  major  types  of  ground  cover.  Short-range  variation  is  also  evident 
related  to  the  water  bodies,  buildings  and  the  more  local  changes  in  ground  cover  and 
drainage.  These  were  described  in  the  previous  final  report. 


Table  2:  Model  parameters  for  the  variograms  computed  for  the  128  by  128  pixel  area 
of  the  SPOT  image 


Variable 

Model  type 

Nugget 

variance 

Sill(l) 

variance 

Range(l) 
pixels  (m) 

Sill(2) 

variance 

Range(2) 

pixels 

(m) 

Nested 

Exponential 

0.0 

227.0 

mvwm 

1 

4.40 

■ 

IB 

Nested 

Spherical 

0.0 

152.2 

mmmm 

II 

91.71 

1BH 

Hermite 

polynomials 

Nested 

Spherical 

0.0 

0.5240 

rnamm 

1 

0.3755 

ISn 

Long-range 

component 

Circular 

0.0 

113.9 

BOH 

IS SI 

Short-range 

component 

Spherical 

0.0 

87.9 

4.30 

(86) 

Low  ffequ 
-ency  1  in  2 

Nested 

Spherical 

0.0 

146.7 

74.82 

■ 

mEEm 

High  ffequ 
-ency  1  in  2 

Pure  nugget 

High  ffequ 
-ency  1  in  4 

Circular 

0.0 

41.7 

2.944 

(58.88) 

High  ffequ 
-ency  1  in  8 

Circular 

0.0 

68.30 

6.690 

(133.8) 
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Odinary  kriged  estimates 


Figure  3:  Pixel  map  of  the  kriged  NER  of  part  of  the  SPOT  image  (128  by  128 
pixels)  for  Fort  A.  P.  Hill 


Factorial  kriging  and  wavelet  analysis  enable  the  different  spatial  scales  to  be  separated 
in  theory.  For  factorial  kriging  this  is  controlled  by  the  variogram  which  describes  the 
variation  present  in  the  data.  For  the  wavelet  analysis  this  is  controlled  by  the 
resolution  level  which  is  controlled  by  the  octave  bands  and  as  a  consequence  is  more 
arbitraiy.  One  aspect  of  future  research  is  to  consider  how  the  variogram  could  be  used 
to  guide  the  factoring  process  that  controls  the  wavelet  multiresolution  analysis. 

Figure  4  is  a  pixel  map  of  the  kriged  estimates  of  the  long-range  component  of  the 
variation  filtered  using  the  variogram.  The  large  scale  variation  is  related  to  the  main 
relief  features.  The  band  of  dark  colours  in  the  North  and  central  part  of  the  map  are 
damper  and  wetter  areas,  and  the  lighter  ones  the  spurs,  upper  slopes  and  built  areas. 
This  map  could  be  used  effectively  to  guide  future  sampling.  If  the  end-user  is 
interested  in  retrieving  this  level  of  information  then  a  suitable  sampling  interval  can  be 
chosen  using  the  range  of  the  variogram.  A  sampling  interval  of  200  m  would  be 
adequate  to  ensure  that  this  resolution  of  variation  is  identified. 

The  pixel  map  of  the  short-range  variation  (Figure  5)  shows  the  detail  that  is  also 
evident  in  Figure  3,  but  less  clearly  so.  The  lakes  are  recovered  well  by  this  resolution. 
The  dark  patches  in  the  bottom  left  hand  comer  (1  to  20  on  the  x-axis  and  55  to  60  on 
the  y-axis),  in  the  central  area  (64  to  90  on  the  x-axis  and  115  to  125  on  the  y-axis), 
and  at  the  top  of  the  map  (45  to  70  on  the  x-axis  and  158  to  180  on  the  y-axis).  The 
road  running  N-S  is  also  evident  extending  N  along  longitude  100  (on  this  map).  The 
other  short-range  structures  probably  relate  to  changes  in  local  drainage  conditions  and 
vegetation.  For  many  surveys  recovering  this  intensity  of  variation  at  a  scale  of  about 
120  m  would  require  too  much  sampling.  A  sampling  interval  of  50  m  to  60  m  would 
be  needed  to  resolve  this  short-range  variation.  If  a  sampling  scheme  of  about  200  m 
were  recommended  in  relation  to  the  long-range  variation  this  information  on  short- 
range  variation  would  be  lost.  These  maps  enable  us  to  demonstrate  to  the  end-user  the 
extent  of  information  that  is  likely  to  be  lost  by  adopting  the  coarser  sampling. 
Sampling  between  60  m  and  200  m  would  be  of  little  benefit  because  most  of  the 
short-range  variation  would  not  be  identified  and  sampling  at  less  than  200  m  would  be 
inefficient  to  identify  the  long-range  variation.  Variograms  were  computed  from  the 
estimates  of  the  long-range  (Figure  6a)  and  the  short-range  (Figure  6b)  components. 
They  recover  the  spatial  scale  of  the  variation  quite  well,  but  both  variograms  were 
difficult  to  model  satisfactorily. 

For  the  first  wavelet  analysis  the  level  of  resolution  was  26.  The  coefficients  were 
derived  as  described  earlier.  The  low  frequency  and  the  high  frequency  coefficients 
were  reconstructed  by  the  inverse  wavelet  transform,  which  restored  each  of  the  64  by 
64  sets  coefficients  to  the  size  of  the  original  data  set.  These  are  shown  as  pixel  maps 
and  should  be  compared  with  the  appropriate  kriged  and  the  filtered  maps.  Figures  3  to 
5.  In  addition  variograms  were  computed  for  each  of  these  reconstructions:  low 
frequency  (Figure  7a)  and  high  frequency  (Figure  7b). 

The  low  frequency  reconstruction,  Figure  8,  is  very  similar  to  the  ordinary  kriged 
output  for  the  image.  Figure  3.  It  is  important  to  remember  that  the  ordinary  kriged 
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Figure  5:  Pixel  map  of  the  short-range  component  of  the  variation  in  NIR  of 
part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Long-range  component 


Figure  6:  a)  Variogram  of  the  long-range  component,  and  b)  variogram  of  the 
short-range  component  of  the  variation  in  NIR  of  part  of  the  SPOT  image  (128 
by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  7:  a)  Variogram  of  the  low-frequency  component,  and  b)  variogram  of  the 
high-frequency  component  of  the  variation  in  NER  of  part  of  the  SPOT  image 
(128  by  128  pixels)  for  Fort  A.  P.  Hill  from  a  wavelet  analysis  at  a  resolution  of  1 
in2 
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map  was  made  from  estimates  using  all  of  the  data,  whereas  the  low  frequency  wavelet 
reconstruction  used  the  1  in  2  sample,  i.e.  25%  of  the  original  data.  Both  the  long-  and 
short-range  components  of  the  variation  are  evident,  although  there  has  been  some  loss 
of  detail  in  the  short-range  variation.  For  example  the  road  is  less  clear  in  Figure  8  than 
in  Figure  3.  The  variogram  computed  from  the  low  frequency  reconstruction,  Figure 
7a,  was  very  similar  to  the  variogram  of  the  raw  data,  Figure  2a.  Hence  the  spatial 
structure  at  both  scales  has  been  retained  at  this  level  of  resolution.  The  most 
surprising  finding  was  that  related  to  the  high  frequency  reconstruction.  The  map  of 
the  high  frequency  component  (Figure  9)  does  not  appear  to  reflect  the  kind  of 
variation  present  in  the  map  of  the  short-range  component  from  factorial  kriging. 
However,  when  we  examined  them  in  detail  there  is  some  weak  evidence  of  the  lakes, 
which  are  very  clear  in  Figure  5,  and  the  road.  The  variogram  computed  from  these 
data  is  pure  nugget,  Figure  7b,  which  means  that  the  high  frequency  components  are 
noise  at  this  level  of  resolution;  they  contain  no  spatial  structure.  The  latter  is  all 
retained  in  the  low  frequency  reconstruction. 

To  determine  whether  we  could  retrieve  the  long-  and  short-range  components  using 
wavelets  we  explored  the  next  resolution,  25,  in  effect  a  sampling  of  1  in  4  (or  1  pixel 
in  16).  Figure  10  shows  the  low  frequency  reconstruction.  There  is  still  long-  and 
short-range  variation  evident,  although  the  short  range  variation  is  becoming  less 
distinct;  for  example  the  road  and  the  lakes  are  still  visible  but  their  margins  are  less 
clearly  defined.  Figure  1 1  shows  the  pixel  map  for  the  average  of  the  high  frequency 
reconstruction  and  it  is  clear  that  there  is  more  of  the  short-range  component  of  the 
variation  evident.  The  variogram  of  the  high  frequency  reconstruction  now  shows 
some  structure.  Figure  12a.  Table  2  gives  the  model  parameters  of  this  variogram. 

The  low  frequency  reconstruction  of  the  1  in  8  resolution  24  now  shows  the  long-range 
component  of  the  variation  identified  by  factorial  kriging.  Figure  13.  This  resolution  is 
fairly  close  to  the  short-range  component  of  the  variogram,  i.e.  6.5  pixels,  and  this 
level  of  variation  appears  to  have  been  filtered  out  now.  So  it  seems  that  once  the 
resolution  of  the  short-range  structure  has  been  reached  the  effect  was  to  remove  the 
short-range  variation.  The  map.  Figure  14,  of  the  high  frequency  reconstruction  now 
shows  some  of  the  features  evident  in  the  kriged  map  of  the  short-range  component  of 
the  variation.  In  particular  the  lakes  are  evident.  The  variogram  computed  from  the 
average  of  the  three  high  frequency  reconstructions.  Figure  12b,  shows  clear  evidence 
of  structure  and  the  range  of  spatial  correlation  described,  6.69  pixels,  is  close  to  the 
short-range  component  of  the  variogram  of  the  original  data. 


Summary 

It  is  clear  that  factorial  kriging  works  well  with  multiresolution  data.  The  main  reason 
for  this  is  that  the  filtering  is  controlled  by  the  variogram  which  is  a  function  of  the 
data  being  analysed.  It  is  a  valuable  method  for  directing  future  sampling  for  ground 
surveys  because  it  can  show  what  degree  of  variation  is  likely  to  be  recovered.  The 
multiresolution  analysis  using  wavelets  produces  a  different  outcome.  At  the  first 
resolution  the  high  frequency  components  remove  the  noise,  i.e  spatially  uncorrelated 
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Wavelet  reconstruction  for  1  in  2  selection 
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Figure  8:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
at  a  resolution  of  1  in  2 
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Figure  9:  Pixel  map  of  the  average  reconstruction  of  the  high  frequency  wavelets 
from  the  wavelet  analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels) 
for  Fort  A.  P.  Hill  at  a  resolution  of  1  in  2 
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Wavelet  reconstruction  for  the  1  in  4  selection 
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Figure  10:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
at  a  resolution  of  1  in  4 
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Figure  11:  Pixel  map  of  the  average  high  frequency  reconstruction  from  the 
wavelet  analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A. 
P.  Hill  at  a  resolution  of  1  in  4 
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Figure  12:  a)  Variogram  of  the  high-frequency  component  for  the  1  in  4 
resolution,  and  b)  variogram  of  the  high-frequency  component  for  the  1  in  8 
from  the  wavelet  analysis  of  NIT  at  Fort  A.  P.  Hill 
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Wavelet  reconstruction  for  1  in  8  selection 
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Figure  13:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
at  a  resolution  of  1  in  8 
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High  frequency  reconstruction  1  in  8 
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Figure  14:  Pixel  map  of  the  average  high  frequency  reconstruction  from  the 
wavelet  analysis  of  NER  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A. 
P.  Hill  at  a  resolution  of  1  in  8 
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Data  reconstruction 


The  128  by  128  pixels  were  sampled  by  taking  one  pixel  in  every  two  for  each  row  and 
column  (or  a  sample  of  1  in  4)  which  matches  the  first  level  of  resolution  in  the  wavelet 
analysis  26,  one  pixel  in  every  four  for  every  row  and  column  (or  a  sample  of  1  in  16), 
a  wavelet  resolution  of  25,  and  one  pixel  in  every  eight  for  each  row  and  column  (or  a 
sample  of  1  in  64),  a  wavelet  resolution  of  24.  The  low  frequency  wavelet  coefficients 
were  inverted  to  reconstruct  the  image  as  before.  For  this  analysis  pixel  maps  have 
been  prepared  so  that  the  information  can  be  compared  more  directly.  Kriged  estimates 
were  made  to  coincide  with  the  original  data  points  for  each  data  set  using  the 
variogram  model  from  the  frill  set  of  data.  These  maps  are  also  shown  as  pixel  maps. 

To  evaluate  the  accuracy  of  the  estimates  by  the  wavelet  reconstruction  and  kriging 
every  value  was  compared  with  the  original  values  of  NIR.  First  the  differences  were 
calculated  between  the  estimates  and  original  values  for  both  analyses  and  for  the  three 
sub-samples,  and  these  are  shown  as  pixel  maps  (Figures  16  and  17,  21  and  22,  26  and 
27).  The  statistical  distribution  of  these  differences  or  errors  has  also  been  determined 
and  these  are  shown  as  histograms  (Figures  18,  23  and  28).  In  addition  the  mean 
squared  differences  or  mean  squared  error  (mse)  was  calculated  (Table  3). 


Results 

The  results  were  not  entirely  what  we  expected  and  we  have  been  making  sure  that  the 
kriging  program  and  analyses  have  been  correct.  From  the  theory  of  geostatistics  we 
should  expect  that  the  kriged  estimates  would  have  the  smallest  mse,  but  they  do  not 
for  any  of  the  analyses.  It  was  this  that  led  us  to  explore  the  differences  in  more  detail 
to  try  to  gain  insight  into  the  results  from  the  two  methods. 


Sample  of  1  in  2 

The  pixel  maps  for  the  low  frequency  wavelet  reconstruction  and  kriging  from  the  1  in 
2  data.  Figures  8  and  15,  respectively  appear  to  be  very  similar  to  each  other.  The 
slight  ‘spottiness’  evident  on  the  kriged  map  is  because  punctual  kriging  was  used  and 
this  is  a  true  estimator  returning  the  value  at  the  data  points.  Table  3  gives  the  mean 
squared  errors  for  both  methods.  That  for  the  wavelets  is  less.  The  maps  of  errors  or 
comparisons.  Figures  16  (wavelet)  and  17  (kriging),  show  a  similar  pattern  in  general. 
However,  the  differences  between  them  help  to  explain  why  the  mse  is  greater  for  the 
kriged  estimates  than  for  the  wavelet  reconstruction.  There  are  large  differences 
associated  with  the  lakes  where  there  are  clearly  marked  local  changes  associated  with 
boundaries  in  the  variation.  This  is  evidence  of  local  non-stationarity  which  violates  the 
assumptions  of  kriging.  Wavelets  are  known  to  be  suitable  for  dealing  with  local  non- 
stationarity,  and  these  results  support  this.  Kriging  has  the  largest  absolute  differences 
and  there  are  more  of  them  than  for  the  low  frequency  wavelet  reconstruction. 
However,  compared  to  the  number  of  pixels  in  the  data  these  larger  differences  are  few 


26 


Figure  15:  Pixel  map  of  the  kriged  estimates  for  the  1  in  2  sample  of  NIR  of  part 
of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  16:  Pixel  map  of  the  comparisons  between  the  kriged  estimates  for  the  1  in 
2  data  with  the  original  NIR  values  of  part  of  the  SPOT  image  (128  by  128  pixels) 
for  Fort  A.  P.  Hill 
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Comparisons  for  the  wavelet  reconstruction  1  in  2 


Figure  17:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  2  data  with  the  original  NIR  values  of  part  of 
the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 


a) 


Wavelet  1  in  2 


Figure  18:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1 
in  2  sampling  for  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A. 
P.  Hill 
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compared  with  the  many  much  smaller  differences  for  the  majority  of  the  estimates. 

To  explore  the  reasons  for  the  results  in  more  detail  the  histograms  of  the  differences 
were  examined.  Figure  18  a  and  b  are  the  histograms  of  the  differences  for  the  wavelet 
and  kriging  analyses,  respectively.  It  is  clear  that  punctual  kriging,  which  is  a  true 
estimator  at  the  data  points  has  a  larger  number  of  small  errors  than  the  wavelet 
analysis.  However,  this  is  not  consistent  as  the  number  of  data  points  retained  is 
reduced. 


Sample  of  1  in  4 

Figure  19  shows  the  result  of  kriging  this  sub-sample.  It  is  evident  that  much  of  the 
short-range  variation  has  been  lost  even  though  the  variogram  of  the  full  data  set  was 
used.  This  map  is  similar  to  that  for  the  long-range  component.  The  map  of  the  low 
frequency  reconstruction.  Figure  20,  shows  more  of  the  short-range  variation  and 
appears  to  be  much  more  accurate  visually  than  the  kriged  map.  The  maps  of  the 
differences,  Figures  21  and  22,  appear  to  be  similar  overall,  but  closer  examination 
shows  that  the  patches  where  the  differences  are  large  for  kriging  are  more  extensive 
than  those  for  the  wavelet  analysis.  The  values  of  the  MSEs  for  kriging  and  the 
wavelet  analysis  in  Table  3  also  suggest  that  kriging  performs  worse  than  wavelets. 
The  histograms,  Figure  23  a  (kriging)  and  b  (wavelets),  suggest  that  more  of  the  kriged 
values  have  smaller  differences  from  the  original  values  than  those  for  the  wavelet 
reconstruction.  However,  the  number  of  large  errors  is  also  greater  for  the  kriged 
values. 


Sample  of  1  in  8 

Figure  24  shows  the  result  of  kriging  this  sub-sample.  It  is  evident  that  much  more  of 
the  detail  in  the  variation  has  been  lost.  The  pattern  that  is  returned  is  coarse  and  no 
longer  reflects  even  the  long-range  component  of  the  variation  as  accurately.  Figure 
25  for  the  low  frequency  reconstruction  from  the  wavelet  analysis  also  shows  how  the 
detail  has  been  lost.  The  maps  of  the  differences,  Figures  26  and  27  are  again  similar, 
but  as  before  where  the  differences  are  greatest  for  the  kriged  differences  (Figure  25) 
so  their  extent  is  also  more  extensive.  It  seems  from  the  MSEs  Table  3  that  as  the  data 
become  more  sparse  and  separated  by  greater  distances  that  kriging  loses  power  in 
comparison  with  the  wavelet  analysis.  The  histograms,  Figure  28  a  (kriging)  and  b 
(wavelets),  suggest  that  the  wavelet  analysis  has  performed  better  with  this  sub-set  of 
the  data.  In  the  central  part  of  the  distribution  there  is  little  difference  between  the 
differences  for  wavelet  analysis  and  kriging,  but  there  seem  to  be  many  more  large 
errors  for  kriging  than  the  wavelets. 
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Kriged  estimates  for  1  in  4  selection 


Figure  19:  Pixel  map  of  the  kriged  estimates  for  the  1  in  4  sample  of  NIR  of  part 
of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Wavelet  reconstruction  for  the  1  in  4  selection 
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Figure  20:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
at  a  resolution  of  1  in  4 
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Comparisons  for  the  kriged  estimates  for  1  in  4 
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Figure  21:  Pixel  map  of  the  comparisons  between  the  kriged  estimates  for  the  1  in 
4  data  with  the  original  NER  values  of  part  of  the  SPOT  image  (128  by  128  pixels) 
for  Fort  A.  P.  Hill 
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Comparisons  for  the  wavelet  reconstruction  1  in  4 


Figure  22:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  4  data  with  the  original  NIR  values  of  part  of 
the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Wavelet  1  in  4 


Figure  23:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1 
in  4  sampling  for  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A. 
P.  Hill 
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Figure  25:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
at  a  resolution  of  1  in  8 
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Figure  26.  Pixel  map  of  the  comparisons  between  the  kriged  estimates  for  the  1  in 
4  data  with  the  original  NIR  values  of  part  of  the  SPOT  image  (128  by  128  pixels) 
for  Fort  A.  P.  Hill 
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Comparisons  for  the  wavelet  reconstruction  1  in  8 
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Figure  27:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  8  data  with  the  original  NIR  values  of  part  of 
the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  28:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1 
in  8  sampling  for  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A. 
P.  Hill 
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Summary 


The  histograms  of  the  differences  are  perhaps  the  most  illuminating  part  of  this 
analysis.  It  seems  that  we  need  to  explore  more,  but  that  kriging  performs  well  when 
fewer  data  have  been  removed  than  the  wavelet  analysis.  It  also  suggest  that  the  end- 
user  can  be  provided  with  some  insight  to  enable  them  to  choose  which  is  appropriate 
for  their  needs.  It  seems  that  for  the  1  in  2  and  lin  4  data  sub-sets  more  of  the  errors 
are  small  for  kriging  than  for  wavelets,  but  that  the  overall  error  is  least  for  the  wavelet 
analysis.  The  latter  is  clearly  more  successful  at  retaining  the  transition  features  present 
which  kriging  will  not  do  well.  Again  what  does  the  end  user  want? 

Another  thing  that  seems  to  emerge  from  this  investigation  is  that  the  variogram  could 
be  used  to  choose  an  optimal  subset  of  the  data,  based  on  the  distance  between  the 
values.  With  the  1  in  2  sample  both  the  long-range  and  short-range  components  of  the 
variation  were  restored  as  we  should  expect  from  the  correlation  structures  in  the 
variogram:  the  distance  between  the  pixels  was  less  than  the  range  of  the  short-range 
component.  With  the  1  in  4  sample  only  the  long-range  structure  is  successfully 
restored.  If  that  is  what  is  required  then  this  can  be  chosen  in  a  way  that  is  driven  by 
the  data  using  the  variogram. 

It  is  interesting  to  note  that  the  means  of  the  kriged  reconstructed  values,  Table  3,  are 
close  in  each  case  to  the  mean  of  the  original  data,  Table  1.  The  variances  for  the 
kriged  values  decrease  as  the  sampling  intensity  decreases  and  is  evidence  of  the 
smoothing  of  the  variation  that  occurs  with  kriging.  However,  the  variance  of  the 
kriged  values  for  the  1  in  2  sample  is  closer  to  the  original  variance  than  any  of  the 
other  analyses.  The  wavelet  analysis  retains  the  variance  better  as  the  sampling 
intensity  decreases.  Geostatistical  simulation  would  probably  perform  even  better  in 
terms  of  retaining  the  variance  in  the  data  and  this  method  should  also  be  compared 
with  wavelet  analysis  in  the  fixture. 
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Abstract 

Wavelet  analysis  operates  locally  and  can  describe  a  wide  range  of  frequencies 
simultaneously  and  filter  them  by  multi-resolution  analysis.  Kriging  analysis  also  filters 
spatial  variation  at  different  resolutions.  We  compare  the  effectiveness  of  wavelets  and 
factorial  kriging  for  exploring  nested  variation  in  a  SPOT  image.  In  addition  both 
wavelets  and  kriging  can  be  used  to  restore  image  data  after  compression.  We  compare 
the  reliability  of  the  restorations  from  the  two  approaches. 

The  near  infrared  (NIR)  waveband  of  part  of  a  SPOT  image  covering  Fort  A.  P.  Hill 
in  Virginia  was  used  for  these  analyses.  The  region  is  on  the  dissected  Piedmont  area 
of  the  eastern  United  States.  An  area  of  128  by  128  pixels  was  selected  from  the  scene 
for  analysis.  The  experimental  variogram  was  computed  and  modelled  by  a  nested 
spherical  function  with  correlation  structures  of  about  6.5  pixels  and  21  pixels.  The 
variogram  and  factorial  kriging  separated  the  two  main  spatial  features  present.  The 
low-frequency  component  from  the  wavelet  analysis  contained  the  spatial  structure. 
The  long-range  component  became  evident  as  the  resolution  decreased.  The  high- 
frequency  components  removed  only  the  uncorrelated  variation  and  we  could  not 
retrieve  the  short-range  component. 

The  image  was  sampled  so  that  one  in  every  four  pixels  was  retained,  one  in  every  16 
and  one  in  every  64.  Using  the  variogram  model  for  the  full  set  of  data  values  were 
estimated  at  the  former  data  points  by  ordinary  kriging.  The  low-frequency  wavelet 
transform  for  these  resolutions  was  inverted  so  that  the  missing  values  were  restored. 
The  restored  values  from  both  analyses  were  compared  with  the  original  values  and  the 
mean  squared  differences  (MSD)  computed.  For  all  resolutions  the  MSD  was  smaller 
for  the  wavelet  reconstruction.  However,  the  MSD  proved  somewhat  misleading  when 
frequency  distributions  of  the  errors  were  compared.  They  suggested  that  wavelets  are 
more  able  to  deal  with  the  local  fluctuations  present  in  the  image  and  with  local  non- 
stationarity  than  kriging,  but  that  for  the  majority  of  points  the  kriged  estimates  have  a 
smaller  error. 

The  paper  will  be  illustrated  with  maps  of  the  results,  and  we  shall  suggest 
improvements  for  restoring  images  by  kriging. 


